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ANNUAL  INTERIM  SCIENTIFIC  TECHNICAL  REPORT 


Grant : 

AF0SR-83-0150,  Algebra  Institute,  University  of  California,  Santa  Barbara, 
Santa  Barbara,  California  93106 

Title: 

Stability  analysis  of  finite  difference  schemes  for  hyperbolic  systems,  and 
problems  in  applied  and  computational  linear  algebra 

Period: 

May  1,  1983  -  April  30,  1984 

Principal  Investigators:  Moshe  Goldberg,  Marvin  Marcus 

(Numbered  items  refer  to  format  specified  by  AFOSR  for  Annual  Technical 
Reports) 

1.  Summary: 

The  research  is  concerned  with  two  principal  related  areas.  Stability 
analysis  of  finite  difference  schemes  for  hyperbolic  initial-boundary  value  prob¬ 
lems  lead  to  an  investigation  of  bounds  for  matrix  norms  and  condition 
numbers.  The  aim  of  this  research  is  to  provide  a  better  understanding  of  tools 
used  in  the  numerical  analysis  of  such  hyperbolic  systems.  Approximation 
schemes  lead  to  systems  of  linear  algebraic  equations  and  the  stability  of  such 
schemes  depends  on  the  eigenvalues  and  singular  values  of  the  associated 
matrices.  Thus,  in  certain  aspects,  the  analysis  of  a  finite  difference  scheme  is  a 
problem  in  numerical  linear  algebra.  Eigenvalue  localization  and  inequalities  for 
matrix  norms  are  two  pertinent  areas  of  classical  research  in  this  field.  The 
investigators  have  used  methods  from  convex  analysis,  the  theory  of  inequali¬ 
ties,  classical  linear  and  multi-linear  algebra  and  numerical  range  theory  in 
attacking  these  problems.  This  project  should  contribute  to  better  understand¬ 
ing  of  advanced  computational  techniques,  and  to  the  improvement  of  basic 
mathematical  tools  often  used  in  numerical  analysis  and  other  fields  of  applied 
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mathematics. 
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2,3.  Research  objectives  and  status  of  research 

The  research  completed  by  Moshe  Goldberg  under  Air  Force  Grant  AFOSR- 
83-0150  during  May  1983  —  April  1984,  consists  of  the  following  two  topics: 

1.  Convenient  Stability  Criteria  for  Difference  Schemes  of  Hyperbolic  Initial- 
Boundary  Value  Problems 

Consider  the  first  order  system  of  hyperbolic  partial  differential  equations 

du(x,t)/dt  =  Adu(z,t)/dx  +  Bu(z,t)  +  f(x,t),  *i0,  t  fe  0, 

where  u(x,t)  is  the  unknown  vector;  A  a  hermitian  matrix  of  the  form 
A  =  Ai  ©  Ag,  where  A\  is  negative  definite  and  Az  is  positive  definite;  and  / (x,t ) 
is  a  given  vector.  The  problem  is  well  posed  in  £2(0,  °°)  if  initial  values 

u(x,t)~  u°(x)  e  X.2(0,®),  x  3s  0, 
and  boundary  conditions 

«i(0.0  =  Suz(0.t)  +g(t),  13  0. 

are  prescribed.  Here,  U]  and  uz  are  the  inflow  and  outflow  parts  of  u 
corresponding  to  the  partition  of  A,  and  5  is  a  coupling  matrix. 

In  the  past  year,  E.  Tadmor  and  M.  Goldberg,  [16],  have  succeeded  in 
obtaining  new,  easily  checkable  stability  criteria  for  a  wide  class  of  finite 
difference  approximations  for  the  above  initial-boundary  value  problem.  The 
difference  approximations  consist  of  a  general  difference  scheme  —  explicit  or 
implicit,  dissipative  or  not,  two-level  or  multi-level  -  and  boundary  conditions  of 
a  rather  general  type. 

Attention  is  restricted  to  the  case  where  the  outflow  boundary  conditions 
are  translatory,  i.e.,  determined  at  all  boundary  points  by  the  same  coefficients. 
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This,  however,  is  not  a  severe  limitation  since  such  boundary  conditions  are 
commonly  used  in  practice.  In  particular,  when  the  numerical  boundary  con¬ 
sists  of  a  single  point,  the  boundary  conditions  are  translatory  by  definition. 

Throughout  the  paper  [18]  it  is  assumed  that  the  basic  scheme  is  stable  for 
the  pure  Cauchy  problem,  and  that  the  assumptions  which  guarantee  the  vali¬ 
dity  of  the  stability  theory  of  Gustafsson,  Kreiss  and  Sundstrom  [18]  hold.  With 
this  in  mind  the  question  of  stability  for  the  entire  difference  approximation  is 
raised. 

The  first  step  in  the  stability  analysis  was  to  prove  that  the  approximation  is 
stable  if  and  only  if  the  scalar  outflow  components  of  its  principal  parts  are 
stable.  This  reduces  the  global  stability  question  to  that  of  a  scalar  homogene¬ 
ous  outflow  problem  of  the  form 

du/  dt  =  adu/ dx,  a  -  constant  >0.  x  &  0,  f  a:  0 
u(x.O)  =  u°(x),  x  a  0;  u(0,f)  =  0,  t  &  0. 

The  stability  criteria  obtained  in  [18]  for  the  reduced  problem  depend  both 
on  the  basic  difference  scheme  and  on  the  boundary  conditions,  but  very  little 
on  the  interaction  between  the  two.  Such  criteria  eliminate  the  need  to  analyze 
the  intricate  and  often  complicated  interaction  between  the  basic  scheme  and 
the  boundary  conditions,  hence  providing  in  many  cases  convenient  alternatives 
to  the  well  known  stability  criteria  of  Kreiss  [22],  and  of  Gustafsson,  Kreiss  and 
Sundstrom  [18].  It  should  be  pointed  out  that  the  old  scheme-independent  sta¬ 
bility  criteria  in  [13,14]  easily  follow  from  the  present  criteria  in  [18]. 

Having  the  new  criteria  in  [18],  all  the  examples  in  the  previous  papers 
[13,14]  were  reestablished.  For  instance,  if  the  basic  scheme  is  arbitrary  (dissi¬ 
pative  or  not)  and  the  boundary  conditions  are  generated  by  either  the  explicit 
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or  implicit  right-sided  Euler  schemes,  then  overall  stability  is  assured.  For  dis¬ 
sipative  basic  schemes  stability  is  proved  if  the  boundary  conditions  are  deter¬ 
mined  by  either  oblique  extrapolation,  the  Box-scheme,  or  by  the  right-sided 
weighted  Euler  scheme.  These  and  other  examples  incorporate  most  of  the 
cases  discussed  in  recent  literature  [2],  [3],  [13],  [14],  [18],  [19],  [21],  [23],  [27], 
[30-35],  [37]. 

Some  new  examples  appear  in  [16]  as  well.  Among  these  it  is  found  that  if 
the  basic  scheme  is  arbitrary  and  two-level,  then  horizontal  extrapolation  at  the 
boundary  maintains  overall  stability.  Other  stable  cases  occur  when  the  basic 
scheme  is  given  by  either  the  backward  (implicit)  Euler  scheme  or  by  the 
Crank-Nicolson  scheme,  and  the  boundary  conditions  are  determined  by  oblique 
extrapolation.  Such  examples,  where  neither  the  basic  scheme  nor  the  boun¬ 
dary  conditions  are  necessarily  dissipative,  could  not  have  been  handled  by  the 
previous  results  in  [13,14]. 

An  extended  version  of  [16],  which  includes  additional  examples  and 
remarks,  is  now  in  final  stages  of  preparation,  [17]. 

Such  contributions  should  be  helpful  to  applied  mathematicians  and 
engineers  in  better  understanding  and  exploiting  old  and  new  finite  difference 
approximations  to  hyperbolic  systems. 

2.  Submultiplic  ativity  and  Other  Properties  of  Ip  Norms  for  Matrices 

In  the  past  three  years,  E.  G.  Straus  (now  deceased)  and  M.  Goldberg,  [9,10], 
investigated  submultiplicativity  properties  of  norms  and  seminorms  on  operator 
algebras  —  an  important  subject  in  many  fields  of  numerical  analysis  and  applied 
mathematics.  In  this  work  an  arbitrary  normed  vector  space  V  over  the  com¬ 
plex  field  C,  with  an  algebra  L(V)  are  studied.  If  N  is  positive  definite,  i.e., 
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N(A)  >  0  for  all  A  *  0,  then  N  is  called  a  generalized  operator  norm.  If  in  addi¬ 
tion,  N  is  (sub-)  multiplicative,  namely  N(AB)  s  N(A)N(B)  for  all  A,  B  e  L(V), 
then  N  is  called  an  operator  norm  on  L{  V). 

Given  a  seminorm  N  on  1(7)  and  a  fixed  constant  p>  0,  then  obviously 
Np  a  fiN  is  a  seminorm  too.  Similarly,  is  a  generalized  operator  norm  if  and 
only  if  N  is.  In  both  cases,  may  or  may  not  be  multiplicative.  If  it  is,  then  p 
is  said  to  be  a  multiplicativity  factor  for  N. 

Having  these  definitions  the  following  is  proved  in  [9]: 

(i)  If  N  is  a  nontrivial  seminorm  or  a  generalized  operator  norm  on  L{V),  then 
N  has  multiplicativity  factors  if  and  only  if 

p*  =  svp\N(AB)  :  N(A)  =  N(B)  =  lj  <  -. 

(ii)  If  hn  <  then  p,  is  a  multiplicativity  factor  for  N  if  and  only  if  p  at  p*. 
Special  attention  was  given  to  the  finite  dimensional  case  where  it  suffices, 

of  course,  to  consider  Qxn*  the  algebra  of  nxn  complex  matrices.  Following 
Ostrowski,  [28],  in  this  case  the  terms  generalized  matrix  norm  and  matrix 
norm  are  adopted  instead  of  generalized  operator  norm  and  operator  norm, 
respectively.  In  this  case  it  is  proved  that  while  nontrivial,  indefinite  seminorms 
on  never  have  multiplicativity  factors,  generalized  matrix  norms  always 

have  such  factors.  In  the  infinite  dimensional  case,  however,  the  situation  was 
less  decisive,  i.e.,  there  exist  nontrivial  indefinite  seminorms  and  generalized 
operator  norms  on  L(V)  which  may  or  may  not  have  multiplicativity  factors. 

In  both  the  finite  and  infinite-dimensional  cases  it  is  proved  that  if  M  and  N 
are  seminorms  on  L(V)  such  that  M  is  multiplicative,  and  if  rj  f  >  0  are  con¬ 
stants  satisfying 


<U(A)  <  N(A)  <  r)U{A)  for  all  A  e  L(V), 


then  any  n  with  fi  >  i)/  (*  is  a  multiplicativity  factor  for  N. 

Using  these  results  it  is  proved,  for  example,  that  if  V  is  an  arbitrary  Hil¬ 
bert  space  and 

r(A)  =  supf|(Ax,z)|  :  x  e  V,  I* (  =  1{.  A  e  L(V), 

is  the  classical  numerical  radius,  then  ycr  is  an  operator  norm  if  and  only  if 
4.  This  assertion  is  of  interest  since  the  numerical  radius  r  is  perhaps  the 
best  known  nonmultiplicative  generalized  operator  norm  [1,4,15,20,29],  and  it 
plays  an  important  role  in  stability  analysis  of  finite  difference  schemes  for 
multi-space-dimensional  hyperbolic  initial-value  problems  [15,24,25,36]. 

Straus  and  Goldberg  also  investigated  C-numerical  radii  which  constitute  a 
generalization  of  the  classical  numerical  radius  r,  defined  in  [7]  as  follows:  For 
given  matrices  A,  C  G  C^,  the  C-numerical  radius  of  A  is 

rc(A)  =  max]  |  tr(CU*AU) |  :  U  nxn  unitary]. 

In  [7]  (compare  [26]),  it  is  shown  that  r  is  a  norm  on  Cnxn  '*  and  so  has  multipli¬ 
cativity  factors  —  if  and  only  if  C  is  not  a  scalar  matrix  and  tr  C  *  0.  Multiplica¬ 
tivity  factors  for  the  above  rc  were  found  in  [7-10,12]. 

In  the  most  recent  effort,  Straus  and  Goldberg,  [ll],  studied  the  well  known 
Ip  norms 

I A  \p  =  I  I*  ]1/p>  ^  =  (®y )  ^  Gixjf  1  ^  P  ^  80 • 

It  was  shown  by  Ostrowski,  [28],  that  these  norms  are  multiplicative  if  and  only  if 
1  asp  *s  2.  For  p  2  it  is  shown  that  is  a  multiplicativity  factor  for  |  A |p  if  and 
only  if  thus,  in  particular,  obtaining  the  useful  result  that 

|  A  |p  is  a  multiplicative  norm  on  Q,xn. 


Continuing  this  effort,  Goldberg  obtained  [5,6]  the  best  possible  constants 
fi(p,q)  and  /x(g.p)  (for  arbitrary  1  Sp.g  s  “)  such  that  relations  of  the  form 

\AB\p  < n(p%q)\A\p\B\q,  \AB\P  ^Kq.p)\A\q\B\p 

hold  whenever  the  matrix  product  AB  exists.  This  leads  to  the  best  possible  con¬ 
stant  X(p,g)  for  which 

IU4|(p  «A(p,g)M|,.  A  e  Crnyn,  (l^p.g  °°). 

where 

||i4|jp  =  maxj  \Ax  \p  :  x  e  C”,  |x \p  =  lj 

is  the  ordinary  Ip  operator  norm  of  an  arbitrary  mxn  matrix  A.  Such  inequali¬ 
ties  could  be  useful  in  determining  the  power  boundedness  of  a  matrix  --  a  basic 
question  in  stability  analysis. 

The  research  of  Marvin  Marcus  under  Air  Force  Grant  AF0SR-83-0150  during 
May  1983  —  April  1984  is  described  below. 

In  the  mathematical  modelling  of  physical  phenomena,  a  standard  tech¬ 
nique  for  solving  the  resulting  partial  diffferential  equation  boundary  value  prob¬ 
lem  is  to  approximate  the  differential  system  at  a  discrete  set  of  points  by  the 
solution  of  a  linear  system.  In  general,  such  a  system  is  large,  non-symmetric 
and  sparse.  The  Tchebychev  iteration  based  on  Tchebychev  polynomials  can  be 
used  to  solve  non-symmetric  linear  systems  whose  eigenvalues  lie  in  the  right 
half-plane.  Moreover,  many  factorizations  and  splitting  techniques  applied  to 
symmetric  systems  yield  non-symmetric  systems  with  spectra  in  the  right-half 
plane.  Manteuflel  [Numer.  Math.  28,  307-327,  1977]  showed  that  the  Tchebychev 
iteration  for  an  N-square  real  linear  system 
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whose  eigenvalues  lie  in  the  right-half  plane,  can  be  carried  out  with  two  param¬ 
eters,  c,d  which  arise  in  evaluating  Fn(«f/c)  where  Tn(z)  =  cosh(ncosA_1(z ))  is 
the  ntH  Tchebychev  polynomial. 

Manteuflel  proved  that  if  the  convex  hull,  H(A),  of  the  spectrum  of  A  is 
known  (for  normal  matrices  this  is  the  numerical  range  of  A  ),  then  the  parame¬ 
ters  c  ,d  in  the  above  iteration  can  be  chosen  to  be  optimal  in  a  minimax  sense. 
In  a  later  paper  [Numer.  Math.  31,  183-208,  1978]  he  discusses  a  method  of 
estimating  H{A)  during  iteration  from  the  sequence  of  residual  vectors.  He  also 
shows  that  a  power  method  variant  of  the  Tchebychev  procedure  yields  eigen¬ 
value  estimates  that  lie  in  the  numerical  range  F(A).  B.  A.  Carrfe,  L  A. 
Hagemann,  R.  B.  Kellog,  and  others,  have  also  done  work  on  dynamic  estimation 
of  the  optimal  SOR  parameter. 

At  the  30th  Anniversary  Meeting  of  the  Society  for  Industrial  Sc  Applied 
Mathematics  in  1982,  (supported  in  part  by  a  grant  from  the  Air  Force  Office  of 
Scientific  Research),  Dr.  Martin  Schultz  of  the  Yale  University  Department  of 
Computer  Science  posed  the  following  general  question  stemming  from  the  work 
of  Manteuffel  and  its  continuation  by  the  Yale  group. 

Obtain  computational  estimates  of  the  set  H(A)  and  the  field  of  values  (i.e., 
numerical  range)  F(A). 

In  work  currently  underway  by  M.  Marcus  and  M.  Sandy,  computer  codes 
have  been  developed  for  obtaining  explicit  numerical  plots  of  H(A)cF(A)  that 
are  useful  in  determining  the  optimal  c,d  described  above.  Work  is  also  in  pro¬ 
gress  to  obtain  computable  bounds  on  the  discrepancy  between  F(A)  and  H(A). 
This  discrepancy  problem  has  been  considered  earlier  in  various  forms  by  many 
authors:  H.  Wielandt,  B.  N.  Moyls,  I.  Filippenko,  B.  Shure,  M.  Sandy,  K.  Fan,  C.  A. 
Berger,  M.  Newman,  R.  C.  Thompson,  P.  R.  Halmos. 
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Another  of  the  areas  of  research  under  this  grant  is  concerned  with  the 


condition  number  ic(A)  of  a  matrix  A,  defined  by 

x(A)  =  i|A|| 

which  plays  an  important  role  in  error  estimates  and  termination  criteria  for 
numerical  iterative  methods  for  the  solution  of  linear  systems.  More 
specifically,  the  following  problem  was  examined  in  publication  §2  listed  under 
Publications  below. 

How  does  the  condition  number  ic(A-B)  of  the  Hadamard  product 

A-B  =  [ay  by  3 

of  two  matrices  A  =  [ciy]  and  B  =  [by]  .  depend  on  the  norms,  singular 
values,  etc.  of  the  factors  A  and  B1 

In  particular,  work  appearing  in  the  preceding  publication  investigates  this 
and  related  problems  for  the  class  of  von  Neumann  norms.  Recall  that  a  von 
Neumann  norm,  !1A'!|,  also  called  a  unitarily  invariant  norm,  satisfies 

\\uxn- 11*11 

for  all  unitary  matrices  U  and  V.  von  Neumann’s  theorem  states  that  if 
Oi  >  •  ■  •  a  an  are  the  singular  values  of  X,  then  ||A"||  can  be  written  as 


11*11  =  p(a  I . «») 

where  p  is  a  symmetric  gauge  function.  Note  that  the  usual  Ip-  norms  arise 
from  gauge  functions: 


<p(x  i . Xn)  =  max 

«sn 


i/p 
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where  k  is  a  fixed  integer,  1  ^  k  s  n,  and  p  ^  1,  e.g.,  if  k  =  1  and  p  =  l  then  <p  spe¬ 
cializes  to  the  Hilbert  norm  or  maximum  singular  value  of  X 


In  the  preceding  formula  ||x||  is  the  usual  Euclidean  norm.  In  other  words, 


. *n)  =  max  |  xa[\)  I 


which  obviously  satisfies  the  definition  of  a  symmetric  gauge  function. 

It  is  worthwhile  to  note  that  there  exists  a  simple  connection  between  the 
condition  number  k(A)  and  the  Hadamard  product.  It  is  proved  in  publication 
#2  listed  below,  that  for  the  Hilbert  norm  the  following  inequality  is  available: 


\\A-B\\*\\A\\-m. 


In  case  B  -  A then 


\\A-A-'\\*\\A\\-\\A-'\\. 


Thus  for  the  usual  matrix  norm  subordinate  to  the  Euclidean  vector  norm 


k(A)  \\A-A- 


It  follows  easily  that 


k(A)  ^  L 


where  L  is  the  largest  product 


L  =  max  |  ay  by  | , 

w 


and  B  =  A  Hence,  in  particular,  if  a  lower  bound  for  the  modulus  of  a  single 
element  of  A~l  is  available,  say  |6*^0|  as  £o.  then 


« 


k(A)>  l“v0|£o- 


These  preliminary  results  suggest  a  number  of  possible  areas  of  investiga¬ 
tion  which  are  currently  underway; 

Investigate  the  inequality 

\\A-B\\*\\AW\m 

for  the  class  of  von  Neumann  norms  and  obtain  easily  computable  lower 
bounds  on  the  corresponding  condition  numbers. 

For  the  linear  system  Ax  =  6  these  considerations  show  that  the  estimated  rela¬ 
tive  error  in  x  can  be  as  large  as 

(a  +  0). 

where  a  and  0  are  the  relative  errors  in  A  and  b  respectively. 
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